Problem Set 1

alt text

Problem Set 2

Write an R function to factorize a square matrix A into LU or LDU, whichever you prefer. You don’t have to worry about permuting rows of A and you can assume that A is less than 5x5, if you need to hard-code any variables in your code. If you doing the entire assignment in R, then please submit only one markdown document for both the problems.

#First, the easy example from Wk2 Notes:
M <- matrix(c(2,6,1,8),2,2)

elim <- diag(length(M[1,]))

piv.col<-1
row.order<-1

elim[piv.col+row.order,piv.col] <- (-1)*(M[piv.col,piv.col]/M[piv.col+row.order,piv.col])**(-1)

print(elim)
##      [,1] [,2]
## [1,]    1    0
## [2,]   -3    1
#the product of: (1) the inverse of the elimination matrix (2) the elimination matrix and 
#(3) the original matrix should equal the original matrix M:

solve(elim)%*%(elim%*%M) == M
##      [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUE
#Next, the 3x3 example
A <- matrix(c(1,2,3,1,1,1,2,0,1),nrow=3)

#Function that calculates the elimination matrices
elim.matrix<-function(M,piv.col, row.order){
  elim <- diag(length(M[1,]))
  
  elim[piv.col+row.order,piv.col] <- (-1)*(M[piv.col,piv.col]/M[piv.col+row.order,piv.col])**(-1)
  
  return(elim)}

#change original matrix subscript "a, b,..." for each step that a calculation uses the matrix A
A.a <- A

E21 <- elim.matrix(A.a,1,1)
E31 <- elim.matrix(A.a,1,2)

A.b <- E31 %*% E21 %*% A

E32<-elim.matrix(A.b,2,1)

U <- E32 %*% E31 %*% E21 %*% A

L <- solve(E21) %*% solve(E31) %*% solve(E32)

Strategy: build out elimination algorithm that:
1. Defines piv.col (the column of the privot), and row.order, the “row order” (my made up term for whether you’re eliminating the first row below the pivot row, second row, etc..) agorithmically depending on the size of the matrix, not hard-coded as in the above example.
2. Calculates the elimination matrices.
3. Keeps track of the elimination matrices so they may be used to find L and U.

Architecture note: group steps into separate functions. Use apply instead of loops if practical.

A <- matrix(c(1,2,3,1,1,1,2,0,1),nrow=3)

elim.matrix<-function(M,piv.col, row.order){
  elim <- diag(length(M[1,]))
  
  elim[piv.col+row.order,piv.col] <- (-1)*(M[piv.col,piv.col]/M[piv.col+row.order,piv.col])**(-1)
  
  return(elim)}

Strategy for piv.col and row.order: Notice in example there’s a pattern of input that corresponds to the row reduction algorithm’s steps. It can be described as a 2x3 matrix:

steps.3 <- matrix(c(1,1,2,1,2,1),3,2)

#I found the matrices that describe the steps for 4x4 and 5x5 matrices by hand. 
#Since the problem states that the largest matrix tested is 5x5, these could be hard coded (and that seems simplest at this point). 
#However, these matrices do follow a factorial pattern and could be found for any reasonably sized square matrix. 

#for a 4x4
steps.4 <- matrix(c(1,1,1,2,2,3,1,2,3,1,2,1),6,2)

#for a 5x5
steps.5 <- matrix(c(1,1,1,1,2,2,2,3,3,4,1,2,3,4,1,2,3,1,2,1), 10,2)

Rewrite the elim.matrix function to take a row above matrices as arguments.

A <- matrix(c(1,2,3,1,1,1,2,0,1),nrow=3)

elim.matrix<-function(M, step){
  piv.col<-step[1]
  row.order<-step[2]
  elim <- diag(nrow = nrow(M))
  elim[piv.col+row.order,piv.col] <- (-1)*(M[piv.col,piv.col]/M[piv.col+row.order,piv.col])**(-1)
  
  elim <- as.matrix(elim)
  return(elim)}

Now, step through the appropriate “steps” matrix to calculate each elimination matrix:

for(i in (1:nrow(steps.3))){print(is.matrix(elim.matrix(A,steps.3[i,])))}
## [1] TRUE
## [1] TRUE
## [1] TRUE

We need these as seperate objects for more calculations:

elim.list <-vector("list", nrow(steps.3))

for(i in (1:nrow(steps.3))){
  elim.list[[i]]<-elim.matrix(A,steps.3[i,])}

Finally, calculate U and L from the list of elimination matrices. This shows that the LU decomposition is working correctly.

(L <- solve(elim.list[[1]]) %*% solve(elim.list[[2]]) %*% solve(elim.list[[3]]))
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    2    1    0
## [3,]    3    1    1
(U <- elim.list[[3]] %*% elim.list[[2]] %*% elim.list[[1]] %*% A)
##      [,1] [,2] [,3]
## [1,]    1    1    2
## [2,]    0   -1   -4
## [3,]    0   -1   -1
L %*% U == A
##      [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE

Now to test on 4x4 and 5x5 matrices. I did test on random sample coefficients from 0 to 50 and ran into some problems with floating-point precision - an interesting note and someting understand in applications.

#It's important to check that the first coefficient isn't zero since there isn't any swapping of rows. Also, make sure they're invertable...remarkably, when I used set.seed(4092834) it came up with two identical rows in the 4x4 matrix. 

set.seed(224)
B <- sample(0:10, 16, replace=T)
B <- matrix(data = B, 4, 4)
solve(B)
##            [,1]       [,2]       [,3]        [,4]
## [1,]  1.0588235  0.9117647 -3.3823529  0.05882353
## [2,]  0.3823529  0.1764706 -0.7352941 -0.11764706
## [3,]  0.4411765  0.5882353 -1.6176471 -0.05882353
## [4,] -2.3823529 -2.1764706  7.7352941  0.11764706
C <- sample(1:10, 25, replace=T)
C <- matrix(C, 5, 5)
solve(C)
##             [,1]        [,2]         [,3]        [,4]        [,5]
## [1,] -0.37352351  0.36304881  0.239803878 -0.31669267  0.28125696
## [2,]  0.12993091 -0.05170493  0.102741253  0.13104524 -0.24938712
## [3,] -0.13215957 -0.10006686  0.082460441  0.06396256  0.10786717
## [4,]  0.07978605  0.03342991 -0.230220637  0.01872075  0.06641409
## [5,]  0.21060842 -0.15756630 -0.001560062  0.07176287 -0.12636505

4x4

elim.list.4 <-vector("list", nrow(steps.4))

for(i in (1:nrow(steps.4))){
  elim.list.4[[i]]<-elim.matrix(B,steps.4[i,])}

#Hard-coding the final dot products is a bad solution. Need to get this into an apply function...but it is proving difficult to work with a list of matrices. 
(L.4 <- solve(elim.list.4[[1]]) %*% solve(elim.list.4[[2]]) %*% solve(elim.list.4[[3]]) %*% solve(elim.list.4[[4]]) %*% solve(elim.list.4[[5]]) %*% solve(elim.list.4[[6]]) )
##      [,1] [,2] [,3] [,4]
## [1,] 1.00    0 0.00    0
## [2,] 0.75    1 0.00    0
## [3,] 0.50    1 1.00    0
## [4,] 1.25    0 0.75    1
(U.4 <- elim.list.4[[6]] %*% elim.list.4[[5]] %*% elim.list.4[[4]] %*% elim.list.4[[3]] %*% elim.list.4[[2]] %*% elim.list.4[[1]] %*% B)
##      [,1]     [,2]  [,3]    [,4]
## [1,]    8   7.0000  4.00  5.0000
## [2,]    0  -2.2500  7.00  1.2500
## [3,]    0   1.7500 -5.00 -0.7500
## [4,]    0 -10.0625  1.75 -0.6875
L.4 %*% U.4 == B
##      [,1] [,2] [,3] [,4]
## [1,] TRUE TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE TRUE
## [4,] TRUE TRUE TRUE TRUE

5x5

elim.list.5 <-vector("list", nrow(steps.5))

for(i in (1:nrow(steps.5))){
  elim.list.5[[i]]<-elim.matrix(C,steps.5[i,])}


(L.5 <- solve(elim.list.5[[1]]) %*% solve(elim.list.5[[2]]) %*% solve(elim.list.5[[3]]) %*% solve(elim.list.5[[4]]) %*% solve(elim.list.5[[5]]) %*% solve(elim.list.5[[6]]) %*% solve(elim.list.5[[7]]) %*% solve(elim.list.5[[8]]) %*% solve(elim.list.5[[9]]) %*% solve(elim.list.5[[10]]))
##      [,1]      [,2] [,3] [,4] [,5]
## [1,]  1.0 0.0000000 0.00  0.0    0
## [2,]  0.8 1.0000000 0.00  0.0    0
## [3,]  0.8 0.5714286 1.00  0.0    0
## [4,]  0.6 1.2857143 2.00  1.0    0
## [5,]  1.0 0.4285714 1.75  0.9    1
(U.5 <- elim.list.5[[10]] %*% elim.list.5[[9]] %*% elim.list.5[[8]] %*% elim.list.5[[7]] %*% elim.list.5[[6]] %*% elim.list.5[[5]] %*% elim.list.5[[4]] %*% elim.list.5[[3]] %*% elim.list.5[[2]] %*% elim.list.5[[1]] %*% C)
##      [,1]      [,2]      [,3]      [,4] [,5]
## [1,]    5  4.000000  3.000000  8.000000 10.0
## [2,]    0  3.800000 -0.400000  1.600000 -7.0
## [3,]    0 -1.371429  1.828571 -4.314286  2.0
## [4,]    0  4.457143  3.057143 11.771429  0.0
## [5,]    0 -4.240000 -1.780000 -2.730000 -2.5
L.5 %*% U.5 == C
##      [,1]  [,2]  [,3]  [,4]  [,5]
## [1,] TRUE  TRUE  TRUE  TRUE  TRUE
## [2,] TRUE  TRUE  TRUE  TRUE  TRUE
## [3,] TRUE  TRUE  TRUE FALSE FALSE
## [4,] TRUE FALSE FALSE FALSE FALSE
## [5,] TRUE FALSE FALSE FALSE FALSE
#while the exact answer isn't reached, it appears due to floating point errors.
C - L.5 %*% U.5
##      [,1]         [,2]         [,3]         [,4]         [,5]
## [1,]    0 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [2,]    0 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [3,]    0 0.000000e+00 0.000000e+00 8.881784e-16 8.881784e-16
## [4,]    0 1.776357e-15 1.776357e-15 1.776357e-15 1.776357e-15
## [5,]    0 2.664535e-15 3.552714e-15 1.776357e-15 1.776357e-15